map<-st_as_sf(sovereignty110)
par(mar=c(0,0,0,0), mfrow=c(1,1))
plot(map$geometry,col="lightyellow")

dim(map)

sel<-map[map$adm0_a3 %in% c("FRA", "NCL","ATA","ATF","USA", "PRI","CHN","TWN","MAR", "SAH","CHN","TWN","ISR","PSX"),c("sovereignt","sov_a3","type","admin", "adm0_a3","name","note_adm0","iso_a3","wb_a3")]
kable(sel)

map<-st_as_sf(sovereignty110)
par(mar=c(0,0,0,0))
plot(map$geometry,col="lightyellow")

dim(map)

par(mfrow=c(1,2))
plot(map[map$sov_a3 %in% c("ISR","JOR","SYR","LBN","EGY"),]$geometry, col=c("gray80","orange","gray80","gray80","gray80"))
title("Limits of Israël",cex=0.5)
plot(map[map$sov_a3 %in% c("MAR","SAH","DZA","MRT"),]$geometry, col=c("gray70","orange","gray70","lightyellow"))
title("Limits of Morocco")

map<-st_as_sf(countries110)
small<-st_as_sf(tiny_countries110)
par(mar=c(0,0,0,0), mfrow=c(1,1))
plot(map$geometry,col="lightyellow")
plot(small$geometry,col="red", add=T)

map<-st_as_sf(countries50)

par(mar=c(0,0,0,0))
plot(map$geometry,col="lightyellow")

# Load map
map<-st_as_sf(countries110)
map<-map[c("adm0_a3","name")]
names(map)<-c("iso3c","name","geometry")

# Add polygons center
coo<-st_coordinates(st_centroid(map,of_largest_polygon = T))
map$Lon<-coo[,1]
map$Lat<-coo[,2]


# adjust some codes
map$iso3c[map$iso3c=="KOS"]<-"XKX"    # Kosovo
map$iso3c[map$iso3c=="PSX"]<-"PSE"    # Palestinian territories
map$iso3c[map$iso3c=="SDS"]<-"SSD"    # South Sudan

# Save
st_write(map,"data/world_map_low.shp",delete_dsn=T)
saveRDS(map,"data/world_map_low.Rdata")

# Load map
map<-st_as_sf(countries50)
map<-map[c("adm0_a3","name")]
names(map)<-c("iso3c","name","geometry")

# Add polygons center
coo<-st_coordinates(st_centroid(map,of_largest_polygon = T))
map$Lon<-coo[,1]
map$Lat<-coo[,2]


# adjust some codes
map$iso3c[map$iso3c=="KOS"]<-"XKX"    # Kosovo
map$iso3c[map$iso3c=="PSX"]<-"PSE"    # Palestinian territories
map$iso3c[map$iso3c=="SDS"]<-"SSD"    # South Sudan

# Save
st_write(map,"data/world_map_high.shp",delete_dsn=T)
saveRDS(map,"data/world_map_high.Rdata")

str(cat$countries)

pays<-cat$countries[cat$countries$income_level!="Aggregates",c("iso3c", "country","capital_city","longitude","latitude", "region","income_level")]
kable(head(pays))

indic<-cat$indicators
kable(head(indic,3))

world_data <- wb_data(indicator = c("AG.SRF.TOTL.K2","AG.LND.ARBL.HA","SP.POP.TOTL","SP.URB.TOTL","NY.GDP.MKTP.CD", "EN.ATM.CO2E.KT"),
                 return_wide = TRUE,
                 start_date = 1990,
                 end_date = 2019,
                 country ="countries_only")

world_data<-world_data[,-c(1,3)]

# recodage (attention : ordre alphabetique)
names(world_data)<-c("iso3c","date","ARB","SRF","CO2","GDP","POP","URB")

pays<-cat$countries[cat$countries$income_level!="Aggregates",c("iso3c", "country","capital_city","longitude","latitude", "region","income_level")]
tab<-merge(pays,world_data, by="iso3c",all.x=F,all.y=T)
kable(head(tab))

saveRDS(tab,"data/wb_don_1990_2019.Rdata")

# Extract meta
meta<-cat$indicators[cat$indicators$indicator_id %in% c("AG.SRF.TOTL.K2","AG.LND.ARBL.HA","SP.POP.TOTL","SP.URB.TOTL","NY.GDP.MKTP.CD", "EN.ATM.CO2E.KT"),]

# Select information and add personal code
meta<-data.frame(meta[,c(1,2,4,5)])
meta$shortcode<-c("ARB","SRF","CO2","GDP","POP","URB")
meta<-meta[,c(5,1,2,3,4)]

# Display
kable(meta)

# Save meta
saveRDS(meta,"data/wb_don_1990_2019_meta.Rdata")
write.table(meta,"data/wb_don_1990_2019_meta.csv",sep=";",dec=",", row.names = F,fileEncoding = "UTF-8")

# Load data
don<-readRDS("data/wb_don_1990_2019.Rdata")

# Correction terres arables
length(don$ARB[don$date==2016])
don$ARB[don$date==2016]<-don$ARB[don$date==2015]
don$ARB[don$date==2017]<-don$ARB[don$date==2015]
don$ARB[don$date==2018]<-don$ARB[don$date==2015]
don$ARB[don$date==2019]<-don$ARB[don$date==2015]

## Correction superfice
don$SRF[don$date==2019]<-don$SRF[don$date==2018]

# Correction CO2
don$CO2[don$date==2016]<-don$CO2[don$date==2015]
don$CO2[don$date==2017]<-don$CO2[don$date==2015]
don$CO2[don$date==2018]<-don$CO2[don$date==2015]
don$CO2[don$date==2019]<-don$CO2[don$date==2015]

# Sudan and south Sudan area
don$SRF[don$iso3c=="SDN"]<-1731610
don$SRF[don$iso3c=="SSD"]<-619700


# Sudan and south Sudan arable area
don$ARB[don$iso3c=="SDN"]<-19823160
don$ARB[don$iso3c=="SSD"]<-619700*19823160/1731610

# GDP Corée du Nord
don$GDP[don$iso3c=="PRK"]<-don$CO2[don$iso3c=="PRK"]*don$GDP[don$iso3c=="KOR"]/don$CO2[don$iso3c=="KOR"]

# URB Kosovo
don$URB[don$iso3c=="XKX"]<-don$POP[don$iso3c=="XKX"]*don$URB[don$iso3c=="SRB"]/don$POP[don$iso3c=="SRB"]

tabmis<-don %>% select(iso3c,country,date,ARB,SRF,CO2,GDP,POP,URB)  %>%
  gather("ARB", "SRF","CO2","GDP","POP","URB",key = VAR, value = DON) %>%
  mutate(Missing = is.na(DON))

kable(head(tabmis))             

res <- tabmis %>% group_by(VAR, date) %>% 
                  summarise(nbmis=sum(Missing), nb=n(), pct = 100*nbmis/nb) 

p<-ggplot(res) + 
  aes(x=date,y=pct,color=VAR) +
  geom_line() +
  scale_y_continuous("% de pays à valeurs manquantes", breaks = c(0,10,20,30,40,50,60,70,80,90,100))+
  scale_x_continuous("Année", breaks = c(1990, 1995, 2000, 2005, 2010, 2015, 2020))
p

res <- tabmis %>% group_by(iso3c, country) %>% 
                  summarise(nbmis=sum(Missing), nb=n(), pct = 100*nbmis/nb) %>%
                  arrange(-pct)

kable(head(res,20))

# load a data.frame with columns space, time, indic 
estim<-function(df   = don)               # dataframe
 
{


# Select variable
sel<-dcast(df,space~time)
M<-as.matrix(sel[,-1])
rownames(M)<-sel$space

####################################################
########   Estimation Procedure ####################
####################################################

# check for lines or column  filled with missing values and eliminate
check1<-apply(!is.na(M),FUN="sum",1)
check2<-apply(!is.na(M),FUN="sum",2)
M<-M[check1>0,check2>0]



dim(M)


# define dimensions
nr<-nrow(M)
nr
nc<-ncol(M)
nc


# create time matrix with NA
Mt<-matrix(rep(1:nc,nr),nrow=nr,nc=nc,byrow=T)
rownames(Mt)<-rownames(M)
colnames(Mt)<-colnames(M)
Mt[is.na(M)]<-NA
Mt[1:4,1:4]



# identify previous and next available valuer
Mt_min<-Mt
Mt_max<-Mt
for ( i in 1:nr) { 
  for ( j in 1:nc) { 
    if (is.na(Mt[i,j])) {
      Mt_min[i,j]<-max(Mt[i,1:j],na.rm=T)
      Mt_max[i,j]<-min(Mt[i,j:nc],na.rm=T)    
    }  
  }
}
Mt_min[1:4,1:4]
Mt_max[1:4,1:4]


# choose estimation method
M_met<-matrix("OK",nrow=nr,ncol=nc)
rownames(M_met)<-rownames(M)
colnames(M_met)<-colnames(M)
for ( i in 1:nr) { 
  for ( j in 1:nc) { 
    if (is.na(Mt[i,j])) {
      M_met[i,j]<-"IN"
      if (is.infinite(Mt_max[i,j])) {M_met[i,j]<-"EX"} 
      if (is.infinite(Mt_min[i,j])) {M_met[i,j]<-"RE"}  
    }   
  }
}
M_met[1:4,1:4]

###### Estimation of missing values  ####

M_est<-M
str(M_est)

#### step 1: Interpolation ######

for ( i in 1:nr) { 
  for ( j in 1:nc) { 
    if ((M_met[i,j]=="IN")) {
      t0<-Mt_min[i,j]
      t1<-Mt_max[i,j] 
      tacm<-(M[i,t1]/M[i,t0])**(1/(t1-t0))
      M_est[i,j]<-M[i,t0]*(tacm**(j-t0))
    }   
  }
}
M_est[1:4,1:4]


#### step 2: Extrapolation ######

for ( i in 1:nr) { 
  for ( j in 2:nc) { 
    if ((M_met[i,j]=="EX")) {
      N<-M_est
      t0<-j-1
      t1<-j
      N<-M_est[is.na(M_est[,t0])==F,]
      N<-N[is.na(N[,t1])==F,]
      tacm<-(sum(N[,t1])/sum(N[,t0]))**(1/(t1-t0)) 
      M_est[i,j]<-M_est[i,t0]*(tacm**(t1-t0))
    }   
  }
}
M_est[1:4,1:4]


#### step 3: Retropolation ######

for ( i in 1:nr) { 
  for ( j in (nc-1):1) { 
    if ((M_met[i,j]=="RE")) {
      N<-M_est
      t0<-j
      t1<-j+1
      N<-M_est[is.na(M_est[,t0])==F,]
      N<-N[is.na(N[,t1])==F,]
      tacm<-(sum(N[,t1])/sum(N[,t0]))**(1/(t1-t0)) 
      M_est[i,j]<-M_est[i,t1]/(tacm**(t1-t0))
    }   
  }
}
M_est[1:4,1:4]



#######################################################
############# EXPORT RESULTS  #########################
#######################################################

x<-reshape2::melt(M_est)
names(x)<-c("space","time","estim")
y<-reshape2::melt(M_met)
names(y)<-c("space","time","method")
z<-merge(x,y, by=c("space","time"))

return(z)
}

# Estim ARB
df<-data.table(space=don$iso3c,time=don$date,index=don$ARB)
est<-estim(df)[,1:3]
names(est)<-c("iso3c","date","ARB_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)

# Estim SRF
df<-data.table(space=don$iso3c,time=don$date,index=don$SRF)
est<-estim(df)[,1:3]
names(est)<-c("iso3c","date","SRF_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)

# Estim CO2
df<-data.table(space=don$iso3c,time=don$date,index=don$CO2)
est<-estim(df)[,1:3]
names(est)<-c("iso3c","date","CO2_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)

# Estim GDP
df<-data.table(space=don$iso3c,time=don$date,index=don$GDP)
est<-estim(df)[,1:3]
names(est)<-c("iso3c","date","GDP_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)

# Estim POP
df<-data.table(space=don$iso3c,time=don$date,index=don$POP)
est<-estim(df)[,1:3]
names(est)<-c("iso3c","date","POP_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)

# Estim URB
df<-data.table(space=don$iso3c,time=don$date,index=don$URB)
est<-estim(df)[,1:3]
names(est)<-c("iso3c","date","URB_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)

# Save
saveRDS(don,"data/wb_don_1990_2019.Rdata")

sel <- dev %>% filter(date %in% c(1990,2019)) %>% filter(POP > 5000000)
p<-ggplot(sel,aes(x=DEV_ECO, y=..density..))+
          geom_histogram(col="black",fill="lightyellow",bins=15) +
          geom_density(kernel="gaussian",col="red")+
          scale_x_log10("PIB (en $/hab)") +
          scale_y_continuous("pays > 5 M. hab") +
  facet_grid(cols = vars(date))
  
p

sel <- dev %>% filter(date %in% c(1990,2019)) %>% filter(POP > 5000000)
p<-ggplot(sel,aes(x=DEV_URB, y=..density..))+
          geom_histogram(col="black",fill="lightyellow",bins=15) +
          geom_density(kernel="gaussian",col="red")+
          scale_x_continuous("Population urbaine (en %)") +
          scale_y_continuous("pays > 5 M. hab") +
  facet_grid(cols = vars(date))
  
p

sel <- dev %>% filter(date %in% c(1990,2019))  %>% filter(POP > 5000000)
p<-ggplot(sel,aes(x=DEV_DUR, y=..density..))+
          geom_histogram(col="black",fill="lightyellow",bins=15) +
          geom_density(kernel="gaussian",col="red")+
          scale_x_log10("Emissions de CO2 (en t./hab)") +
          scale_y_continuous("pays > 5 M. hab") +
  facet_grid(cols = vars(date))
  
p

sel <- dev %>% filter(date %in% c(1995,2015)) %>%
               filter(POP > 5000000) %>% 
               pivot_wider(id_cols = iso3c, 
                           names_from= date,
                           values_from = DEV_DUR )
names(sel)<-c("Code","Xt1","Xt2")
sel$var<-100*(((sel$Xt2/sel$Xt1)^(1/20))-1)

p<-ggplot(sel,aes(x=Xt1,y=var,label=Code)) + 
        geom_point(size=1,col="red",alpha=0.2) +
        geom_text(size=2,col="black") +
  geom_smooth(method = "lm")
p